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Abstract 


Many Polyakov loop models can be written in a dual formulation which 
is free of sign problem even when a non-vanishing baryon chemical potential 
is introduced in the action. Here, results of numerical simulations of a dual 
representation of one such effective Polyakov loop model at finite baryon 
density are presented. We compute various local observables such as energy 
density, baryon density, quark condensate and describe in details the phase 
diagram of the model. The regions of the first order phase transition and 
the crossover, as well as the line of the second order phase transition, are 
established. We also compute several correlation functions of the Polyakov 
loops. 
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1 Introduction 


The properties of strongly interacting matter at finite temperatures and densities 
remain in the focus of intensive theoretical and numerical studies (see Ref. 
for a recent review). The full understanding of these properties is still far from 
satisfactory, especially at finite baryon chemical potential, due to the famous sign 
problem. Many approaches to solve this problem, partially or completely, have 
been designed during the last decades and a certain progress has been achieved 
within such methods as Taylor expansion and reweighting at small baryon chemical 
potential, simulations at imaginary potential, complex Langevin simulations and 
some others (see, e.g., the reviews pH). 

One of the approaches attempting to fully solve the sign problem relies on 
rewriting the original partition function and important observables in terms of 
different, usually integer valued, degrees of freedom such that the resulting Boltz- 
mann weight is positive definite. Conventionally, all such formulations are referred 
to as dual formulations, though sometimes a different name can be used (e.g., flux 
line representation) EI There are several routes to construct such dual theory for 
non-Abelian lattice models with fermions og, While a dual formulation with 
a positive Boltzmann weight has not yet been constructed for full QCD, positive 
formulations (or formulations where sign problem appears to be very soft) are 
already known for few important cases. One such case refers to the strong cou- 
pling limit of QCD, where the SU(N) lattice gauge theory can be mapped onto 
a monomer-dimer and closed baryon loop model o (for a recent development of 
this direction, see Ref. and references therein). Another important case is rep- 
resented by many effective Polyakov loop models which can be derived from the 
full lattice QCD in certain limits. A dual representation with positive Boltzmann 
weight is known for some SU(N) and U(N) Polyakov loop models [6] (for 
recent advances, see Ref. IO. Two of these versions have been studied numeri- 
cally in (12\[13)[15}. The emphasis in these simulations was put on establishing the 
phase diagram of the model in the presence of the baryon chemical potential and 
on computing local observables which can be obtained by differentiating the dual 
partition function with respect to some of the parameters entering the action of 
the theory. 

Important class of observables not yet computed in dual formulations are the 
correlations of the Polyakov loops. These correlations can be related to screening 
(electric and magnetic) masses at finite temperatures. Understanding the proper- 
ties of such masses would lead to an essential progress in our comprehension of the 
high temperature QCD phase as a whole. While correlations and related masses 
have been subject of numerous and intensive calculations at zero chemical potential 
(see and references therein), it seems to be an extremely difficult problem to 
compute these masses for the real baryon chemical potential with available simula- 
tion methods. So far correlations and screening masses have been computed only 
at imaginary chemical potential in a7. A closely related and intriguing problem is 
the appearance of a hypothetical oscillating phase at finite density [181120]. Such a 
phase is ultimately connected to the complex spectrum of the theory and requires 
computations of long-distance correlations with real baryon chemical potential. 

Here and in a forthcoming paper we study a somewhat different, but equivalent 
dual form of the effective Polyakov loop model presented in [14]. This form of the 


dual representation has been already used by us in for the computation of 
correlation functions related to the three-quark potential. We believe it is well 
suited to address the problem of screening masses at finite densities, at least in 
the framework of the available positive dual formulations. In the present paper 
we describe the Polyakov loop model we work with, its dual representation and 
several observables. We compute also some local observables and reveal the phase 
structure of the model. We shall also present preliminary results for the Polyakov 
loop correlations. In a companion paper we will give a detailed study of screen- 
ing masses, based on the computation of correlation functions and of the second 
moment correlation length at finite density. This will also allow us to draw some 
conclusions about the existence of an oscillating phase in the model. 

We work on a 3-dimensional hypercubic lattice A = Li. with L the linear 
extension and a unit lattice spacing. The sites of the lattice are denoted by 7 = 
x£ = (x1, 22,23), xi € |0, L — 1], while I = (Z,v) is the lattice link in the direction 
V; e, is the unit vector in the direction v and N; is the lattice size in the temporal 
direction. Periodic boundary conditions are imposed in all directions. Let G be 
the SU(N) group and U(x) an element of G, then dU denotes the (reduced) Haar 
measure on G and TrU the fundamental character of G. 

In this paper we shall study an effective 3-dimensional Polyakov loop model 
which describes a (3 + 1)-dimensional lattice gauge theory with one flavor of stag- 
gered fermions. The general form of the partition function of the model is given 
by 


Za(B,m,u; N,N) =Z = f IT Wæ) [[ 2:00) [[ Balm.) . (1) 


Here, B,(8) is the gauge part of the Boltzmann weight and B,(m, u) is the deter- 
minant for static quarks. There are many forms of B,() and By(m, u) discussed 
in the literature. In what follows we use the weight that can be obtained on an 
anisotropic lattice and in the limit of vanishing spatial gauge coupling Bs after 
explicit integration over all spatial gauge fields (see, for instance, Refs. 
and references therein): 


Bg(8) = exp [8 ReTrU(x)TrU'(z + e,)] . (2) 


For SU(N) the effective coupling constant 8 is related to the temporal coupling 
D, by B = 2Dr(6) with 


pa = (gay) + Cel) = So ert 
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where L,(x) is the modified Bessel function and A, refers to the fundamental rep- 
resentation of SU(N) and is equal to A; = ôu. The Boltzmann weight of static 
staggered fermions can be presented as 


By(m, ul = A(m) det [1 + h+U(x)] det [1 +h Ul(z)] , (4) 
where the determinant is taken over group indices and 
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A(m) = hoN f h4 = het T À h= e Nearcsinh m me TT . (5) 


Below we shall use the dimensionless quantities m = mpn/T and u = lpn /T. Then 
for SU(3) one has A(m) = em. 
The resulting Polyakov loop model we work with takes the form 


Z = f [ [av (2) exp 


x I I A(m) det [1 + hj U(æ)] det [1+ h_U'(z)] . (6) 
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In this model the matrices U(x) play the role of Polyakov loops, the only gauge- 
invariant operators surviving the integration over spatial gauge fields and over 
quarks. The integration in (6) is performed with respect to the Haar measure on 
G. The pure gauge part of the SU(N) model is invariant under global discrete 
transformations U(x) + ZU(x), with Z € Z(N). This is the global Z(N) symme- 
try. The quark contribution violates this symmetry explicitly. Another important 
feature of the Boltzmann weight is that it becomes complex in the presence of a 
chemical potential, as it follows from (6). Therefore, the model cannot be directly 
simulated if u is non-zero. 

In the absence of static quarks the Polyakov loop model exhibits a first order 
phase transition at the critical point 8. = 0.274. The global Z(N) symmetry gets 
spontaneously broken above 6.. At finite density the model defined in Eq. (6) 
and its several variations have been studied both numerically via simulation of 
dual formulations and analytically via mean-field approximation 
and via linked cluster expansion [27]. Mean-field and Monte Carlo study are in 
quantitative agreement for the expectation values of energy density and Polyakov 
loop. This allowed to reveal the phase diagram of the model, at least in some 
regions of the parameters 3, h and u. In this paper we confirm the qualitative 
picture found in previous study and give further details on the behavior of local 
observables, including the baryon density and the quark condensate. Also, first 
results for Polyakov loop correlations will be presented. 

The paper is organized as follows. In Section 2 we formulate our dual rep- 
resentation of the model valid for all SU(N) groups and in all dimensions. We 
present also results of an analytical study of the model based on strong coupling 
expansion and mean-field approximation. The phase diagram of the 3-dimensional 
SU(3) model is studied numerically in details in Section 3, where we discuss also 
simulation results for some local observables as the baryon density and the quark 
condensate. In Section 4 we present preliminary results for the Polyakov loop cor- 
relation functions. The summary and outline for future work is done in Section 5. 


2 Dual formulation of the Polyakov loop model 


In this Section we describe the dual form of the partition function (6). This dual 
representation will be used in the next Sections for numerical simulations of the 
model. All details of the derivation can be found in [14]. In the case of one flavor 
of staggered fermions the partition function (6) can be presented, after an exact 


integration over Polyakov loops, as 


DES Mat ann [fAtm)ra(n(a).pla). (7 
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where l;, i = 1,..., 2d are 2d links attached to a site x and 


Ry(n,p)= 3 A. A. Sntkpetran d(o/1*)d(o + qX/1') KRL. (10) 


q=- k,l=0 oF n+k 


The sum over o runs over all partitions of n + k, and d(o/1™) is the dimen- 
sion of a skew representation defined by a corresponding skew Young diagram, 
o +q = (01 +9,...,0n Lol (for more details we refer the reader to Ref. DÉI 
Equation is valid for all SU(N) groups and in any dimension. Clearly, all fac- 
tors entering the Boltzmann weight of are positive. Hence, this representation 
is suitable for numerical simulations. The Kronecker delta-function in expres- 
sion represents the N-ality constraint on the admissible configurations of the 
integer-valued variables s(/) and r(l). This constraint can be exactly resolved only 
in the pure gauge model when h+ = 0. In this case the dual representation has 
been already tested by us on an example of 2-dimensional SU(3) model, where we 
studied correlation functions and three-quark potential (21). 

In the following Sections we study the dual representation via Monte 
Carlo simulations for the 3-dimensional SU(3) model. In this case the function 
Ry(n, p; h+) takes the form 


R3(n,p) = Q3(n + 1,p) (hy +h? +hyhå CNN) (11) 
+ Qa(n,p) (1 +h} +h? + Ain?) + Q3(n,pt 1) (h- + hå t+ hi h_ t+ MARX) 
+ Q3(n+1,p+1) (hyh + hh?) + Q3(n + 2,p)hyh2 + Qa(n, p+ BAER 


The function Q3(n, p) is the result of the group integration and is given by 


Quin) =>, dQ) dA+ la), (12) 
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where d'Al is the dimension of the permutation group S, in the representation A, 
q=(p—n)/N (when q is not an integer vin, p) = 0). 

Important is the fact that both local observables and long-distance quantities 
can be computed with the help of this dual representation. Explicit expressions for 
the correlation functions of the Polyakov loops will be given in Section 4. Below 
we list some local observables which we are going to compute here and which can 
be obtained by taking suitable derivatives of the partition function (7). 


Magnetization and its conjugate 


M= (EU), == ET å (13) 


Susceptibility 
x = BP (((TrU(2))’) — (TU(2))Å) , (14) 


Energy density 


gac ër SC Tal em, (15) 


Baryon density 


1 OlnZ 1 
De a al = 16). (16) 
e Quark condensate 
1 OlnZ 1 


In the last two equations k(x) and I(x) are the summation variables from (10). 

Before discussing results of numerical simulations we would like to present 
some results obtained by simple analytical methods. These results can serve as 
an additional check of numerical data and lead to a better understanding of the 
whole phase diagram of the model and of the behavior of the different expectation 
values. 


2.1 Strong coupling expansion 


The formulation allows a straightforward expansion in powers of 3. The free 
energy of the 3-dimensional SU(3) model can be written as 


F = 3m + In R3(0,0) + Y BY fr . (18) 


k=1 


The first three coefficients fi are given by 


fi = 3 poi Pio 


3 
fo = A (pi + Poz P20 — 22 på; Pio + 5 poz Pio + 5 Poi P20 + 10 pio Poi p11) 
1 
fs = 3 (1168 pb Pio — 960 på Pio Pus — 480 (po1 Po2 Pio + Por Pro P20) 
+ 20 (pos Pio + Po P30) + 150 p11 (p02 Pio + På P20) 
+ 84 po1 Dio (pia + Poz P20) + 60 (po1 Po P12 + Poi Pio P21) 
+ 30 pir (Pio P12 + Por P21) + 15 pen (Pos Pio + Por P12) 
+ 15 po2 (po1 Dan + P10 P21) + Pos P30 + 3P12 P21) , (19) 


where pu = R3(k,1)/R3(0,0) and the coefficients R3 can be easily calculated from 


Eqs. and (12), resulting in 

y=1+h? + h" + hå + 2h? cosh(3y) , 

) = he"? (14 Rh?) + het (1+h? +h’), 

) = he?" (14h?) +he (1th? +h’) , 

) = he™ (1+ h?)? + Re (14h?) , 
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) = 1 + 2h? + 2hf + hê + 2h? cosh(Bp) , 

) = R(3,0) D + h?)? + 2h? cosh(3y) , 

) = 2h"e"" (1 + A?) + het (24 3h? + 2h") , 

) = 2h"e" (1+ h?) + he (24 3h? +20) . (20) 
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All local observables listed above can be obtained from the expansion (18). The 
strong coupling expansion converges, presumably in the region 6 < GU. u), where 
Be(h, u) is the phase transition or crossover point. One expects that in this region 
numerical data agree reasonably well with strong coupling results. For the purposes 
of this paper it was sufficient to consider only the lowest three orders in the strong 
coupling series. Higher orders can be calculated using the linked cluster expansion 
developed in for a similar Polyakov loop model. 


2.2 Mean-field solution 


Another obvious approach which can be used to get qualitative description of the 
model (6) is mean-field approximation. Within the mean-field method one obtains 
an approximate phase diagram of the model. Also, it allows to calculate various 
local observables, as the free energy density, the baryon density and some others. 
We use here one of the simplest mean-field schemes, applied to a similar Polyakov 
loop model in [24126]. In this scheme the mean-field approximation reduces to the 
following replacement: 


Kä ReTrU(x)TrU' (x + eg) — — SEH (wTrU (x) + uTrU'(x)) , (21) 
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where 


1 Oln Zmelu, w) 1 Oln Zmelu, w) 


Sa = = (TrU# (x)) = 22 
(MUE) = gg TL, u= (TUe) = gg TRE, (a) 
The partition function gets the form 
Z = [Zau w), (23) 
Zme(u,w) = A(m) [vex [48 (wTrU + uTrU")| 
x det [1+A,U] det [1+h_U"] . (24) 


Using the integration methods developed in the mean-field partition func- 
tion is presented as 


Zata = 3" RI ABU" oa (25) 
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r,s=0 


and explicit form of R3(r,s) reads 


1 
Ra(r, s) = Kä Cij Q3(r +i, s +j) + Coo Qa(r + 2,5) + Coz Q3(r,s +2), (26) 


i,j=0 
where the coefficients Ci; are given by 


Coo = 2cosh3m + 2 cosh äu Ou =2coshm , Co =e", Cop = eC", 
Cio = 2e™ oeh ien + 2e?" coshm , Co, = 2e” cosh 2m + 2e~7 cosh m . (27) 


The mean-field equations (22), as well as all local observables, can now be com- 
puted numerically and compared with the strong coupling results and numerical 
data. In the pure gauge case, h = 0, one finds a first order phase transition at the 
critical value fs = 0.2615. This value, as well as the magnetization, matches the 
corresponding results obtained earlier by the mean-field method [p425]. Fig. 
compares our numerical data with the strong coupling expansion and mean-field 
results for some typical values h = 0.6,u = 0.5. More mean-field results and 
comparison with simulations will be given in the next Section. 


3 Phase diagram of the model 


3.1 Lattice setup 


To explore the phase structure of the model, we simulate numerically the partition 
function (7). Important ingredient of such simulations is the way we treat the 
triality constraint in (10). In the pure gauge theory, h = 0, this constraint is 
solved exactly in terms of genuine dual variables and the resulting theory can be 
simulated with the usual Metropolis update (21). Instead, at non-zero values of h, 
the function Ry(n, p) in (7) is explicitly expanded in series (11). In this formulation 
every configuration of link variables s(/) and r(l) has non-zero Boltzmann weight, 
thus allowing us to use again the simple Metropolis update algorithm instead 
of more complicated worm-like algorithms usually adopted to probe dual model 
formulations. The values of the function Q3(n,p) are computed beforehand and 
stored in an array, which is then used in simulations. An extra benefit of the 
absence of the triality condition is the possibility to calculate a correlation function 
as the expectation value of a product of one-site observables instead of a product 
over the path connecting the sources (see Section 4). 

Let us mention that our dual formulation allows to simulate all SU(N) models 
on equal footing: the only difference between different N is encoded in the function 
Qx(n,p) which, as said above, can be computed prior to simulations. Thus, the 
present approach can be easily extended to any values of N. The investigation 
of the large-N limit is interesting per sé since it can shed light on the large-N 


Figure 1: Comparison of the observables obtained from mean-field analysis (solid 
lines), strong coupling expansion up to third order in 8 (dashed lines) and simu- 
lation results (dots) at h = 0.6, u = 0.5 on a 16° lattice. Top left: magnetization 
(blue and red lines denote magnetization and its conjugate, respectively); top right: 
energy density; bottom left: baryon density; bottom right: quark condensate. 


QCD phase diagram at finite density, which could exhibit novel phases, such as 
quarkyonic matter (29131). Moreover, it can reveal interesting connections between 
glueball states in SU(N) gauge models for large N and the constituents of dark 
matter (32]. 

For simulations above h = 0.01, we performed 2 - 10° thermalization updates, 
and then made measurements every 100 whole lattice updates, collecting a statis- 
tic of 10° measurements. To estimate statistical error, a jackknife analysis was 
performed at different blocking over bins with size varying from 10 to 104. 

When we move to smaller values of h, we see that transition probabilities of 
the one-spin Metropolis update become smaller. To counter this, while keeping 
the simulation time reasonable, we used the combination of a multihit update with 
an update skipping procedure: we go through each link variable and attempt Hun 
Metropolis updates, but before each attempt we calculate the a priori probability 


that the update is rejected and generate the number mye; of rejects before the first 
acceptance (from the corresponding geometric distribution); in this way we can 
skip nej updates and therefore save on the number of updates left to perform on 
the current link. Then, if we still need to do some updates, we choose the update 
for the link using the probabilities calculated in presumption that the update 
is accepted. Using this technique allows us to vary nnit for different simulation 
parameters in the range 5 - 2000, while keeping the simulation time constant for 
a fixed acceptance rate. 

Below h = 0.01 we performed 2-10° thermalization updates, with measurements 
taken every 50 whole lattice updates, collecting a statistics of 10° measurements. 
The bins size in the jackknife analysis varied from 10 to 10°. 


3.2 Critical behavior 


A clear indication of the presence of different phases can be seen from the inspection 
of the distribution and the scatter plot of the Monte Carlo equilibrium values for 
the absolute value of the magnetization near a transition value of @: for some 
choice of the parameters h and u, we observe two separate peaks and two separate 
spots, respectively, which is suggestive of a first order transition; for other choices 
we see instead just one single peak and a single spot, respectively, as expected for 
a second order transition or a crossover. 

Another indication comes from the behavior of the magnetization and its sus- 
ceptibility versus 6 near the transition for different values of (h, u). Figs. [2/4] show 
that, when A is kept fixed and p is varied from u = 0 to u = 2, the transition soft- 
ens, the “jump” of the magnetization at the pseudo-critical 9 becoming less steep 
and the peak of the susceptibility less pronounced, which suggests that different 
regimes are being explored. According to the finite-size scaling analysis discussed 
below, the three regimes represented in Figs. [24] correspond to first order, second 
order and crossover, respectively. 
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Figure 2: Magnetization (left) and magnetization susceptibility (right) versus 8 at 
h = 0.01 and u = 0 on a 16? lattice, in the vicinity of a first order phase transition. 
The solid lines represent the mean-field estimates. 
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Figure 3: The same as Fig. [2] for h = 0.01 and u = 0.9635, where a second order 
phase transition occurs. Lines and points in blue correspond to the magnetization, 
while those in red correspond to the conjugate magnetization. 


To characterize in a quantitative way the phase structure of the model in the 
(h, )-plane, we performed a standard finite-size scaling analysis on the peak value 
of the magnetization susceptibility x. Although the simulation algorithm does not 
prevent us from considering arbitrarily large lattice sizes, we limit ourselves to 
the results presented here, obtained by high-statistic simulations on lattices with 
linear sizes L = 10 and 16. Larger lattices were also studied for some choices of 
the couplings to investigate the behavior of the correlation functions. We do not 
quote those results here to keep a uniform treatment of the whole parameter space. 
These results will appear in a separate paper. The reason for using relatively small 
lattice sizes is twofold: on one side, we had to explore a wide region in a three- 
parameter space and simulating many volumes for each point in the parameter 
space would have been too expensive; on the other side, drawing a precision phase 
diagram for this model is beyond the scope of the present work, whose main aim 
is rather to show the effectiveness of the suggested model. We determined y for 
several 3 values in the transition region and fitted them to a Lorentzian, thus 
getting the position of the peak, which gives the pseudocritical coupling fpc, and 
its height. Comparing the dependence of the peak height on the lattice size L with 


the scaling law 
XL( bpc) SC ALY , (28) 


we estimated the critical-index ratio y/v and collected all our determinations, as 


many as 202, in "Table |1| More specifically, for each (h, u) pair we extracted y/v 


by the following formula: 
T = logis (aes) (29) 
H 10 VX Lol 

and assigned to each determination a statistical error calculated by standard error 


propagation. 
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Figure 4: The same as Fig. [3] for h = 0.01 and u = 2, where a crossover transition 
occurs. 
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We can see that, within uncertainties, the values of y/v are spread in a range be- 
tween 3., which implies a first order transition, and zero, which holds for crossover, 
passing through the second order 3-dimensional Ising value, y/v = 1.9638(8) [33]. 
These sparse values of y/v are evidently an artifact of the relatively small lattice 
sizes we could simulate. If we could approach the thermodynamic limit, we would 
see that values of y/v concentrate around the values of 3. (first order), 1.9638(8) 
(second order in the 3-dimensional Ising class) and 0. (crossover). This is expected 
since we know that at u = 0 in the pure gauge limit of QCD or for heavy enough 
quark masses there is a whole region of first order deconfinement transitions in 
the m,a-ms plane (the famous Columbia plot), delimited by a line of second order 
critical points in the 3-dimensional Ising class EI thereafter, for lower quark 
masses, the crossover region is met. In the simulations of our effective Polyakov 
loop model at non-zero density toward the thermodynamic limit we should see the 
continuation of the line of second order critical points to non-zero values of the 
chemical potential. For the lattice volumes considered in our study, we are not 
able to make a clear-cut assignment of each choice of the parameters h and u to 
one of the three transition regions. Using the determination of y/v, we tried any- 
how to make this assignment, extending and modifying the three possible options 
(first order, second order and crossover) as in Table|2| This makes no sense in the 
thermodynamic limit, but can be helpful in the present context. In Table |2| we 
introduced a color code, to help identifying at a glance in which of these regions a 
given parameter pair (h, u) falls. 


y/v color phase 


y/v >3 green first order 
2.50 <y/v <3 | light green | more first order than second order 
1.98 < y/v < 2.50 yellow more second order than first order 


1.94 < y/v < 1.98 red very close to second order 

1.85 < y/v < 1.94 brown more second order than crossover 

0.3 < y/v < 1.85 | magenta | more crossover than second order 
re age Et blue crossover 


Table 2: Color code used to characterize the different phases depending on the 
value of y/v. 


In Fig. [5[ left) each parameter pair (h, u) considered in our simulations is rep- 
resented by a colored dot in the (h, u) plane, according to the color code defined in 
Table Dr allowing us to sketch a tentative phase diagram in the right panel of the 
same figure. A different visualization of the general phase diagram is presented in 
Fig. [6] where in a 3d plot y/v values are reported in correspondence of each (h, u) 
pair. 

Fig. |7| shows the values of Da, obtained for each considered pair (h, u). It 
is interesting to note that the variation of pe along the red line (close to the 
second order phase transition) is much smaller than the variation for all other 
points. Hence, red points lie approximately in one plane pce = const, while all 
points associated with other phases lie on a curved surface. This is also seen 
from Table [3} where we collected the values of pe close to the second order phase 
transition. 
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Figure 5: (Left) Assignment of each parameter pair (h, u) to a transition region 
according to the color code of Table|2} (Right) Estimated phase diagram. 


We have not performed simulations in the absence of external field. To find 
the critical value 8, at h = 0, i.e. for the pure gauge theory, and at u = 0., 
we performed a simple fit of the form 6.(h) = Øy + ah, as suggested in EI We 
obtained following values 8, = 0.274991, a = —0.5568. The value of 8, agrees very 
well with the value quoted in the literature, 8, = 0.274 and reasonably well with 
the mean-field result 6, = 0.2615. 

Another important question concerns the shape of the critical line shown in the 
right panel of Fig. [5]in the heavy-dense limit, h + 0, u + co. From the data we 
have, one cannot make unambiguous conclusions about its behavior. Nevertheless, 
data are well fitted by the function He = —alnh+c—bh?, with a = 0.932, c = —3.1, 
b = 2367. This shows that the line of second order phase transition might persist 
in the heavy-dense limit of QCD. 


3.3 Quark condensate and baryon density 


In this subsection we study the behavior of the quark condensate and the baryon 
density in different phases and compare numerical results with mean-field predic- 
tions. In the static approximation for the quark determinant one cannot observe 
the phenomenon of spontaneous breaking of the chiral symmetry. Indeed, even in 
the strong coupling region, using Eq. (18). one can easily obtain for the quark con- 
densate Q = 0 for all um the massless limit h = 1. The same result in this limit 
can be obtained within mean-field approach and from numerical simulations which 
we performed for various values of 8 and u. Nevertheless, we think it might be in- 
structive to see the behavior of the condensate in the three regimes corresponding 
to first and second order transitions and to crossover. Å more traditional observ- 
able to study in the effective Polyakov loop models is the baryon density. In the 
dual formulation the baryon density B and the quark condensate Q are given by 
Eqs. and (17), respectively. We have computed the right-hand sides of these 
equations both numerically and using the mean-field approach, for the same values 
of h and yu. 
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Figure 6: Values of the critical index ratio y/v for each considered choice of the pair 
(h, u). The color code is as in Table 2] The plane at y/v = 1.9638(8) corresponds 
to the y/v value for a second order phase transition in the 3-dimensional Ising 


class . 


The behavior of baryon density and quark condensate as functions of 6 depends 
strongly on the phase of the system. Left panels of Figs. [8] and [9] show the typical 
behavior of Q and B for h = 0.01 and various values of the chemical potential 
corresponding to different phases of the model. These phases are characterized as 
before by the values of y/v and can be approximately read off from Fig. [6Jor, more 
precisely from the "Table (UL One observes a rapid change of both Q and B at the 
first order phase transition. This rapid change becomes smoother and smoother 
when parameters are gradually changed toward the second order line and then to 
the crossover regime. The right panels of the same Figures compare Monte Carlo 
data with the mean-field predictions in the three regions. One can conclude that 
mean-field reproduces numerical simulations with good accuracy. 

The quark condensate as a function of 3 is shown in Fig. [LO| left) for vanishing 
value of the chemical potential and several values of h. The right panel of the 
same figure demonstrates the approach to the saturation of the baryon density at 
zero quark mass, h = 1. 

As follows from Fig. UL the baryon density B is a decreasing function of 8 at 
sufficiently large h = 0.6, while Fig. [9]shows that for small h values (large mass) B 
is an increasing function of 9. This conclusion is supported by all three methods 
of calculations used in this paper. Therefore, there should exists some value of the 
quark mass where the density changes its qualitative behavior. We have not tried 
to estimate this value. 


4 Large-distance behavior of the correlations 
To see the impact of a non-zero chemical potential on the correlation function 


behavior, we calculated the two-point correlation functions for several values of 
parameters. We considered six kinds of the correlation functions: 
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Figure 7: Values of pc for each considered choice of the pair (h, u). The color 
code is as in Table OU For the sake of readability, red points have been connected 
by a broken line and the plane bpe = 0.266334 was drawn. 


Pan(r) = (TrU(0)TrU(r)) , P(r) = (Re Tr U(0) Re Tr U(r)) , 
hal) = (ru) re), T-i(r) = (Re Tr U(0) Im Tr U(r)) , 
Tarr Oru), (hø (ba Trëtt force Ers (80) 


In the dual formulation the correlation functions can be written as 


Pate) = (0100) Fatale) + LAO) 


Rs(n(0),p(0))_—Raln(r), PC) 
Fu) eat 


(31) 


These formulas work for r > 0. For r = 0 both shifts to the n, p variables happen 
at one point, so only one ratio remains. The correlations [’,,., I; and Ii; can be 
obtained as linear combinations of Dyn, Png and Taa. 

The expressions become unusable when h = 0, and can have a bad con- 
vergence properties for very small h, or very large u values. We have checked by 
comparing the numerical results with the strong coupling expressions for small 6 
values, that the results can be relied on for h > 0.005 and u < 3. 

Since we work at non-zero h, the average traces can become non-zero, intro- 
ducing a constant term into the correlation function even in the disordered phase. 


17 


h H Boe y/v 
0.001 3.36 0.266708(1) | 1.96(2 
0.003 | 2.25 | 0.266724(3) | 1.97(2 


0.01 | 0.963 | 0.26648(1) | 1.95(3 
0.01 | 0.96325 | 0.26649(1) | 1.95(3 
0.01 | 0.9635 | 0.26649(1) | 1.94(3 
0.01 | 0.964 | 0.266473(4) | 1.95(2 
0.01 | 0.96425 | 0.26648(1) | 1.97(3 
0.01 | 0.9645 | 0.26647(1) | 1.98(3 
0.015 0.3 | 0.26625(1) | 1.94(3 
0.0154 | 0.2 | 0.26624(1) | 1.96(4 


( 
( 
( 
( 
( 
( 
( 
( 
0.0155 | 0.08 | 0.266312(5) | 1.95(2 
0.0155 | 0.1 | 0.26630(1) | 1.95( 
0.0156 | 0.0 | 0.266303(4) | 1.94( 
0.0156 | 0.02 | 0.266295(4) | 1.96( 
( 
( 
( 
( 
( 
( 
( 
( 
( 
( 
( 
( 


0.0156 | 0.025 | 0.26629(1) | 1.98(4 
0.0156 | 0.0275 | 0.26631(1) | 1.95(2 
0.0156 | 0.03 | 0.266293(5) | 1.94(3 
0.0157 | 0.0 | 0.26625(1) | 1.97(3 
0.0157 | 0.02 | 0.266247(1) | 1.97(1 
0.0157 | 0.05 | 0.266235(1) | 1.98(3 
0.0158 | 0.0 | 0.266198(3) | 1.95(2 
0.0158 | 0.04 | 0.266191(4) | 1.95(3 
0.01585 | 0.0 | 0.26615(1) | 1.94(4 


0.01585 | 0.01585 | 0.266171(4) | 1.96(3 
0.01585 | 0.02 | 0.26617(1) 
0.01585 | 0.03 | 0.26616(1) 


Table 3: Values of 6pc and y/v for h and u belonging to the region “very close to 
second order”. These values correspond to the bigger red points in Fig. 


Due to that, we introduce the subtracted correlation functions, subtracting the 
corresponding average trace inside the correlations: 


Lulla BM", 
Paal?) = Par = (Tr U) (Tr Ut) ? 
Paasub(T) = Paa(r) — (Tr Ut)? . (32) 


For these subtracted correlations we expect an exponential decay, 


T(r) = Aa l (33) 
at least in the disordered phase. 

Samples of correlation function behavior in different regions of the phase dia- 
gram are shown in Figs. [1 1] [I2Jand[13] One can see that, indeed, both in disordered 
and ordered phases, the correlations decay exponentially. While the mass gap, cor- 
responding to the slope of the plots, remains the same for Han, Pra, Laa, Lu, and 
I,; correlation functions, it is much larger for the T; correlation function. Also the 
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Figure 8: Behavior of the quark condensate versus 6 for fixed h = 0.01 and for 
different values of u on a 16° lattice. Left panel: results of simulations in the 
vicinity of a phase transition. Right panel: comparison of mean-field analysis and 
numerical results in three different phase regimes. The color code of Table |2] is 
used here and below to differentiate the phases of the system. 


mass gap for the I;; correlation remains more or less constant, and in particular 
does not vanish in the vicinity of the phase transition. 

The difference in mass gaps can be explained by noting that at u = 0 I, 
and I’;; correspond to the color-magnetic and color-electric sectors having different 
mass gaps my and mg, My < Mp (see (17). While at non-zero pi these two 
sectors should mix, so all the correlators we study should decay with my, it is 
possible that in our case the mixing is small causing the real large distance mass 
gap for the color-electric sector to be visible only on distances larger than the ones 
for which we have reliable results. 

In the ordered phase we observe an increase of the correlation function slope 
(at the same 8 values) with the increase of u, implying that the mass gap grows 
with u. This means an increase of the screening effects: at finite density non-zero 
u pushes system deeper in the deconfined phase. This is in qualitative accordance 
with the results of obtained with imaginary p. 

We will address these questions in more detail in a future work, which is under 
preparation. For now, we can note that at least in the disordered phase also 
the second moment correlation length for the imaginary-imaginary correlations is 
substantially different from the one for the real-real correlations. 


5 Summary 


Revealing the phase diagram of QCD at finite temperature and non-zero baryon 
chemical potential, as well as the nature of strong interacting matter at high 
temperatures, remains one of the important challenges of high energy physics. 
In spite of enormous efforts, many aspects of these problems are still far from 
unambiguous resolution. The several methods developed to solve the sign problem 
in finite density QCD have their own advantages and drawbacks. Dual formulations 
of effective Polyakov loop models for heavy-dense QCD, as the one used here, 
solve the sign problem completely, but are restricted so far to strong coupled 
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Figure 9: Behavior of the baryon density versus (@ for fixed h = 0.01 and for 
different values of u on a 16 lattice. Left panel: results of simulations in the 
vicinity of a phase transition. Right panel: comparison of mean-field analysis and 
numerical results in three different phase regimes. 


regions in the case of non-Abelian gauge theories. Even in this region one can 
get valuable information about the phase structure of non-Abelian models and 
behavior of various observables. The study presented here is in the same spirit of 
Refs. [12]13][15], though we have used a different dual formulation of the Polyakov 
loop model, built in Ref. (14. To corroborate our findings we have also compared 
in many cases simulation results with the strong coupling expansion of the dual 
model and with the mean-field analysis. Let us briefly recapitulate our main 
results. 


e The phase diagram of the model was studied in great details. We have 
classified three regions in the parameter space of the model (8, h, p) according 
to the type of the critical behavior: first or second order phase transition, 
or crossover. The values of the ratio of critical indices y/v is different in 
different regimes. 


e As main observables we computed expectation values of the Polyakov loop 
and its conjugate, the baryon density and the quark condensate. 


e Our dual formulation allows us to compute correlation functions of the 
Polyakov loops. We have presented some preliminary results for such corre- 
lations at non-zero chemical potential. 


e It is interesting to note that the mean-field results agree very well with 
numerical simulations both at zero and non-zero u. Also, the mean-field 
results are in good qualitative agreement with a similar analysis in Ref. (24). 


The overall qualitative picture of the phase diagram and the behavior of all 
observables fully agree with the picture described in Refs. [13]13]. 

All observables considered in this work have shown sensitivity to the chemical 
potential. The general trend is that when u is increased, they exhibit a less 
steep variation across transition when the coupling 8 (which corresponds to the 
temperature in the underlying QCD theory) is increased. Qualitatively, one can 
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Figure 10: Left panel: behavior of the quark condensate versus £ for fixed u = 0 
and for different values of h on a 16° lattice. Right panel: baryon density for h = 1 
and 6 = 0.1 on a 16° lattice. 


say that increasing u plays effectively the same role as a reduction of the quark 
mass. 

The most important direction for the future work is a detailed study of the 
different correlations of the Polyakov loops and the extraction of screening (electric 
and magnetic) masses at finite chemical potential. Also, an investigation of the 
oscillating phase and the related complex masses can be accomplished within our 
dual formulation. All these problems will be addressed in a companion paper 
which is currently under preparation. 
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Figure 11: Behavior of the correlation functions on a 20° lattice for different values 
of parameter 3 (h = 0.008, u = 0.9635). Top: disordered phase. Middle: near 
the first order phase transition point. Bottom: ordered phase. This legend applies 
also to the next figures. 
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B=0.26649, h= 0.01, u= 0.9635 


Figure 12: The same as Figure|11|for h = 0.01, u = 0.9635 (region of the second 
order phase transition). 
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Figure 13: The same as Figure [11] for h = 0.01, u = 2 (crossover region). 
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